# Replication File for Seljan and Gronke, 2022, "Happy Birthday You Get
# To Vote". 
#
# ANALYSIS: Oregon AVR analysis, OR 2018 election
#
sink("PSRM AVR OR 2018 Replication_log.txt", append=F, type="output")

library(equivtest)
library(data.table)
library(dplyr)
library(lubridate)
library(mosaic)
library(stringr)
library(AER)
library(broom)
library(tidyverse)
library(tidymodels)
library(stargazer)
library(radiant)
library(openxlsx)
options(scipen = 100, digits = 3)


# Variable Creation -------------------------------------------------------

OR_2018<-readRDS(file="Oregon/OR2018_anon.RDS")

##read registration dates as lubridates
OR_2018$birth_year<-OR_2018$BIRTH_DATE.x
OR_2018$lub_regdate<-mdy(OR_2018$EFF_REGN_DATE.x)
OR_2018$reg_year<-year(OR_2018$lub_regdate)
OR_2018$reg_day_of_year<-yday(OR_2018$lub_regdate)

#Dummy variables for general election registration years
OR_2018$regyear2018<-ifelse(OR_2018$reg_year==2018, 1, 0)
OR_2018$regyear2016<-ifelse(OR_2018$reg_year==2016, 1, 0)
OR_2018$regyear2014<-ifelse(OR_2018$reg_year==2014, 1, 0)
OR_2018$regyear2012<-ifelse(OR_2018$reg_year==2012, 1, 0)



# Compute age at November 2018 General Election accounting for leap years
OR_2018<-OR_2018 %>% mutate(age_at_2018election = case_when(BNOV6==1 ~ 2018-birth_year, 
                                                            BNOV6==0 ~ 2018-birth_year-1))

OR_2018<-OR_2018 %>%
  mutate(age1823=ifelse(age_at_2018election>=18 & age_at_2018election<=23, 1, 0),
         age2430=ifelse(age_at_2018election>=24 & age_at_2018election<=30, 1, 0),
         age3140=ifelse(age_at_2018election>=31 & age_at_2018election<=40, 1, 0),
         age4160=ifelse(age_at_2018election>=41 & age_at_2018election<=60, 1, 0),
         age61plus=ifelse(age_at_2018election>=61, 1, 0))

# Limit data to registered voters 18 or over at time of election
OR_2018<-OR_2018 %>%
  filter(age_at_2018election<=110 & age_at_2018election>17)  


# Proxy for new registrations
OR_2018<-OR_2018 %>%
  mutate(new18_v3=ifelse(reg_year==2018 & general16!="NO" & general16!="YES", 1,0))


#Compute dichotomous voting record.
OR_2018$voteyes_2018<-ifelse(OR_2018$general18=="YES",1,0)
OR_2018$voteyes_2016<-ifelse(OR_2018$general16=="YES" ,1,0)
OR_2018$voteyes_2018<-OR_2018$voteyes_2018 %>% replace_na(0)
OR_2018$voteyes_2016<-OR_2018$voteyes_2016 %>% replace_na(0)

#Create Dummies for party status
OR_2018$democrat<-ifelse(OR_2018$PARTY_CODE.x=="DEM",1,0)
OR_2018$republican<-ifelse(OR_2018$PARTY_CODE.x=="REP",1,0)
OR_2018<-OR_2018 %>% mutate(third_party=case_when(PARTY_CODE.x=="REP" ~ 0, 
                                                  PARTY_CODE.x=="DEM" ~ 0,
                                                  PARTY_CODE.x=="NAV" ~ 0,
                                                  TRUE~1))

OR_2018$hasparty<-ifelse(OR_2018$PARTY_CODE.x=="NAV", 0,1)




#Creates Instrumented variable (registered by Oct 16). 
deadline2018<-yday(ymd(20181016))
OR_2018$reg2018_preOct16<-ifelse(OR_2018$reg_day_of_year<=deadline2018 & OR_2018$reg_year==2018,1,0)
OR_2018$reg_preOct16<-ifelse(OR_2018$lub_regdate<=ymd(20181016),1,0)


# Limit data to those with birthdays 42 days prior and 42 days subsequent to AVR Cutoff 
# Coding accounts for leap years. 

OR_2018<-OR_2018 %>%
  mutate(samplewindow = if_else(AAUG14==1 & BNOV6==1, 1, 0) )

#subset of data = registered in 2018 in sample window
sample_OR2018_reg2018<-OR_2018 %>%
  filter(reg_year==2018 & samplewindow==1)

#Create variable equivalent to regyear2018 but named for multi-file use
OR_2018$reg_duringAVR<-OR_2018$regyear2018

reg_duringAVR<-OR_2018 %>%
  filter(reg_year==2018)


# Table 3 -----------------------------------------

table(sample_OR2018_reg2018$BSEP25)

# All registrations (not new)
## 1 = Wald Estimation for subset registered during AVR in birthday window
prop1_1<-prop.test(reg2018_preOct16~BSEP25, data=sample_OR2018_reg2018, success=1)
firststage1<-prop1_1$estimate[2] - prop1_1$estimate[1]
prop2_1<-prop.test(voteyes_2018~BSEP25, data=sample_OR2018_reg2018, success=1)
effect1<-prop2_1$estimate[2] - prop2_1$estimate[1]
wald1<-effect1/firststage1

prop1_1
firststage1
prop2_1
effect1
wald1

sample_OR2018_reg2018 %>%
  group_by(BSEP25) %>%
  filter(!(is.na(reg2018_preOct16))) %>%
  filter(!(is.na(BSEP25))) %>%
  summarize(n=n(),
         proportion=mean(reg2018_preOct16))

iv1_wald<-ivreg(voteyes_2018 ~ reg2018_preOct16 | BSEP25 , data=sample_OR2018_reg2018)
summary(iv1_wald, diagnostics = TRUE)




# 2018 objects for Table 5 -----------------------------------------------------------------


iv1OR18<-ivreg(voteyes_2018 ~ reg2018_preOct16 + populous_county+white+black+hispanic+asian+democrat+republican+third_party+age_at_2018election +female + nogender | BSEP25 +populous_county+white+black+hispanic+asian+democrat+republican+third_party+age_at_2018election + female + nogender, data=sample_OR2018_reg2018)

ivnewOR18<-ivreg(voteyes_2018 ~ reg2018_preOct16 + populous_county+white+black+hispanic+asian+democrat+republican+third_party+age_at_2018election +female + nogender | BSEP25 +populous_county+white+black+hispanic+asian+democrat+republican+third_party+age_at_2018election + female + nogender, data=subset(sample_OR2018_reg2018, new18_v3==1))  

iv1OR18.sum<- summary(iv1OR18, vcov. = function(x) vcovHC(x, type="HC0"), diagnostics=T)
ivnewOR18.sum<- summary(ivnewOR18, vcov. = function(x) vcovHC(x, type="HC0"), diagnostics=T)

iv1OR18.sum
ivnewOR18.sum

iv1OR18.rob<- coeftest(iv1OR18, function(x) vcovHC(x, type="HC0"))
ivnewOR18.rob<- coeftest(ivnewOR18, function(x) vcovHC(x, type="HC0"))

save(iv1OR18, ivnewOR18, iv1OR18.rob, ivnewOR18.rob, ivnewOR18.sum, iv1OR18.sum, file = "Oregon/Object Files 2018/OR 18 iv main robust.RData")


# 2018 - Table 6 -----------------------------------------------------------------

######LIMIT TO NEW
sample_OR2018_reg2018<-sample_OR2018_reg2018 %>%
 filter(new18_v3==1)


iv18_23<-ivreg(voteyes_2018 ~ reg2018_preOct16 + populous_county+white+black+hispanic+asian+democrat+republican+third_party+age_at_2018election +female + nogender | BSEP25 +populous_county+white+black+hispanic+asian+democrat+republican+third_party+age_at_2018election + female + nogender, data=subset(sample_OR2018_reg2018, age_at_2018election>=18 & age_at_2018election<=23))  
iv24_30<-ivreg(voteyes_2018 ~ reg2018_preOct16 + populous_county+white+black+hispanic+asian+democrat+republican+third_party+age_at_2018election +female + nogender | BSEP25 +populous_county+white+black+hispanic+asian+democrat+republican+third_party+age_at_2018election + female + nogender, data=subset(sample_OR2018_reg2018, age_at_2018election>=24 & age_at_2018election<=30))  
iv31_40<-ivreg(voteyes_2018 ~ reg2018_preOct16 + populous_county+white+black+hispanic+asian+democrat+republican+third_party+age_at_2018election +female + nogender | BSEP25 +populous_county+white+black+hispanic+asian+democrat+republican+third_party+age_at_2018election + female + nogender, data=subset(sample_OR2018_reg2018, age_at_2018election>=31 & age_at_2018election<=40))  
iv41_60<-ivreg(voteyes_2018 ~ reg2018_preOct16 + populous_county+white+black+hispanic+asian+democrat+republican+third_party+age_at_2018election +female + nogender | BSEP25 +populous_county+white+black+hispanic+asian+democrat+republican+third_party+age_at_2018election + female + nogender, data=subset(sample_OR2018_reg2018, age_at_2018election>=41 & age_at_2018election<=60)) 
iv61plus<-ivreg(voteyes_2018 ~ reg2018_preOct16 + populous_county+white+black+hispanic+asian+democrat+republican+third_party+age_at_2018election +female + nogender | BSEP25 +populous_county+white+black+hispanic+asian+democrat+republican+third_party+age_at_2018election + female + nogender, data=subset(sample_OR2018_reg2018, age_at_2018election>=61))  

ivfemale<-ivreg(voteyes_2018 ~ reg2018_preOct16 + populous_county+white+black+hispanic+asian+democrat+republican+third_party+age_at_2018election  | BSEP25 +populous_county+white+black+hispanic+asian+democrat+republican+third_party+age_at_2018election , data=subset(sample_OR2018_reg2018, female==1))
ivmale<-ivreg(voteyes_2018 ~ reg2018_preOct16 + populous_county+white+black+hispanic+asian+democrat+republican+third_party+age_at_2018election  | BSEP25 +populous_county+white+black+hispanic+asian+democrat+republican+third_party+age_at_2018election , data=subset(sample_OR2018_reg2018, female==0 & nogender==0)) 
ivnogender<-ivreg(voteyes_2018 ~ reg2018_preOct16 + populous_county+white+black+hispanic+asian+democrat+republican+third_party+age_at_2018election  | BSEP25 +populous_county+white+black+hispanic+asian+democrat+republican+third_party+age_at_2018election , data=subset(sample_OR2018_reg2018, nogender==1)) 

ivpop<-ivreg(voteyes_2018 ~ reg2018_preOct16 + populous_county+white+black+hispanic+asian+democrat+republican+third_party+age_at_2018election +female + nogender | BSEP25 +populous_county+white+black+hispanic+asian+democrat+republican+third_party+age_at_2018election + female + nogender, data=subset(sample_OR2018_reg2018, populous_county==1)  )
ivnonpop<-ivreg(voteyes_2018 ~ reg2018_preOct16 + populous_county+white+black+hispanic+asian+democrat+republican+third_party+age_at_2018election +female + nogender | BSEP25 +populous_county+white+black+hispanic+asian+democrat+republican+third_party+age_at_2018election + female + nogender, data=subset(sample_OR2018_reg2018, populous_county==0))  

ivblack<-ivreg(voteyes_2018 ~ reg2018_preOct16 + populous_county+white+black+hispanic+asian+democrat+republican+third_party+age_at_2018election +female + nogender | BSEP25 +populous_county+white+black+hispanic+asian+democrat+republican+third_party+age_at_2018election + female + nogender, data=subset(sample_OR2018_reg2018, black==1))  
ivwhite<-ivreg(voteyes_2018 ~ reg2018_preOct16 + populous_county+white+black+hispanic+asian+democrat+republican+third_party+age_at_2018election +female + nogender | BSEP25 +populous_county+white+black+hispanic+asian+democrat+republican+third_party+age_at_2018election + female + nogender, data=subset(sample_OR2018_reg2018, white==1))  
ivhispanic<-ivreg(voteyes_2018 ~ reg2018_preOct16 + populous_county+white+black+hispanic+asian+democrat+republican+third_party+age_at_2018election +female + nogender | BSEP25 +populous_county+white+black+hispanic+asian+democrat+republican+third_party+age_at_2018election + female + nogender, data=subset(sample_OR2018_reg2018, hispanic==1))  
ivasian<-ivreg(voteyes_2018 ~ reg2018_preOct16 + populous_county+white+black+hispanic+asian+democrat+republican+third_party+age_at_2018election +female + nogender | BSEP25 +populous_county+white+black+hispanic+asian+democrat+republican+third_party+age_at_2018election + female + nogender, data=subset(sample_OR2018_reg2018, asian==1))  


gaze.coeft <- function(x, col="Std. Error"){
  stopifnot(is.list(x))
  out <- lapply(x, function(y){
    y[ , col]
  })
  return(out)
}

gaze.lines.ivreg.diagn <- function(x, col="statistic", row=1:3, digits=1){
  stopifnot(is.list(x))
  out <- lapply(x, function(y){
    stopifnot(class(y)=="summary.ivreg")
    y$diagnostics[row, col, drop=FALSE]
  })
  out <- as.list(data.frame(t(as.data.frame(out)), check.names = FALSE))
  for(i in 1:length(out)){
    out[[i]] <- c(names(out)[i], round(out[[i]], digits=digits))
  }
  return(out)
}

ivnewOR18.summ <- summary(ivnewOR18, vcov. = function(x) vcovHC(x, type="HC0"), diagnostics=T)
ivnewOR18.rob  <- coeftest(ivnewOR18, function(x) vcovHC(x, type="HC0"))


iv18_23.sum<- summary( iv18_23, vcov. = function(x) vcovHC(x, type="HC0"), diagnostics=T)
iv24_30.sum<- summary( iv24_30, vcov. = function(x) vcovHC(x, type="HC0"), diagnostics=T)
iv31_40.sum<- summary( iv31_40, vcov. = function(x) vcovHC(x, type="HC0"), diagnostics=T)
iv41_60.sum<- summary( iv41_60, vcov. = function(x) vcovHC(x, type="HC0"), diagnostics=T)
iv61plus.sum<- summary( iv61plus, vcov. = function(x) vcovHC(x, type="HC0"), diagnostics=T)
ivfemale.sum<- summary( ivfemale, vcov. = function(x) vcovHC(x, type="HC0"), diagnostics=T)
ivmale.sum<- summary( ivmale, vcov. = function(x) vcovHC(x, type="HC0"), diagnostics=T)
ivnogender.sum<- summary( ivnogender, vcov. = function(x) vcovHC(x, type="HC0"), diagnostics=T)
ivpop.sum<- summary( ivpop, vcov. = function(x) vcovHC(x, type="HC0"), diagnostics=T)
ivnonpop.sum<- summary( ivnonpop, vcov. = function(x) vcovHC(x, type="HC0"), diagnostics=T)
ivwhite.sum<- summary( ivwhite, vcov. = function(x) vcovHC(x, type="HC0"), diagnostics=T)
ivblack.sum<- summary( ivblack, vcov. = function(x) vcovHC(x, type="HC0"), diagnostics=T)
ivhispanic.sum<- summary( ivhispanic, vcov. = function(x) vcovHC(x, type="HC0"), diagnostics=T)
ivasian.sum<- summary( ivasian, vcov. = function(x) vcovHC(x, type="HC0"), diagnostics=T)

iv18_23.rob<- coeftest( iv18_23, function(x) vcovHC(x, type="HC0"))
iv24_30.rob<- coeftest( iv24_30, function(x) vcovHC(x, type="HC0"))
iv31_40.rob<- coeftest( iv31_40, function(x) vcovHC(x, type="HC0"))
iv41_60.rob<- coeftest( iv41_60, function(x) vcovHC(x, type="HC0"))
iv61plus.rob<- coeftest( iv61plus, function(x) vcovHC(x, type="HC0"))
ivfemale.rob<- coeftest( ivfemale, function(x) vcovHC(x, type="HC0"))
ivmale.rob<- coeftest( ivmale, function(x) vcovHC(x, type="HC0"))
ivnogender.rob<- coeftest( ivnogender, function(x) vcovHC(x, type="HC0"))
ivpop.rob<- coeftest( ivpop, function(x) vcovHC(x, type="HC0"))
ivnonpop.rob<- coeftest( ivnonpop, function(x) vcovHC(x, type="HC0"))
ivwhite.rob<- coeftest( ivwhite, function(x) vcovHC(x, type="HC0"))
ivblack.rob<- coeftest( ivblack, function(x) vcovHC(x, type="HC0"))
ivhispanic.rob<- coeftest( ivhispanic, function(x) vcovHC(x, type="HC0"))
ivasian.rob<- coeftest( ivasian, function(x) vcovHC(x, type="HC0"))


stargazer(iv18_23,  iv24_30,  iv31_40,  iv41_60,  iv61plus,  ivfemale,  ivmale, 
          ivpop,  ivnonpop,  ivwhite,  ivblack,  ivhispanic,  ivasian, 
          keep="reg2018_preOct16",
          dep.var.labels =NULL, digits=2, 
          column.labels= c("Aged 18-23", "Aged 24-30", "Aged 31-40", 
                           "Aged 41-60", "Aged 61 Plus", "Female", "Male",  
                           "Populous County", "Nonpopulous County", "White", 
                           "Black", "Hispanic", "Asian"),
          se = gaze.coeft(list(iv18_23.rob,  iv24_30.rob,  iv31_40.rob,  
                               iv41_60.rob,  iv61plus.rob,  ivfemale.rob,  
                               ivmale.rob,  ivpop.rob,  ivnonpop.rob,  
                               ivwhite.rob,  ivblack.rob,  ivhispanic.rob,  
                               ivasian.rob)), df=FALSE, omit.stat=c("adj.rsq","ser"),
          add.lines = gaze.lines.ivreg.diagn(list( iv18_23.sum,  iv24_30.sum,  
                                                   iv31_40.sum,  iv41_60.sum,  
                                                   iv61plus.sum,  ivfemale.sum,  
                                                   ivmale.sum, ivpop.sum,  
                                                   ivnonpop.sum,  ivwhite.sum,  
                                                   ivblack.sum,  ivhispanic.sum,  
                                                   ivasian.sum), row=1), out="Oregon/Object Files 2018/OR18_subgroupiv.tex")


# Note: Table 6 is a combination of results from OR18_subgroupiv.htm, 
# OR16_subgroupiv.htm, and CA18_subgroupiv.htm.


# 2018 objects for Figure 2 ----------------------------------------------------------------


#######Equivalence tests Hildago and Hartmann

dem0<-reg_duringAVR$democrat[reg_duringAVR$BSEP25==0]
dem1<-reg_duringAVR$democrat[reg_duringAVR$BSEP25==1]
dem_equiv<-fisher.binom.two.sided(dem0, dem1, eps_sub=.01, alpha=0.05)
dem_equiv<-c(1, dem_equiv$odds.ratio, dem_equiv$equiv.CI[1], dem_equiv$equiv.CI[2], "Democrat")

rep0<-reg_duringAVR$republican[reg_duringAVR$BSEP25==0]
rep1<-reg_duringAVR$republican[reg_duringAVR$BSEP25==1]
rep_equiv<-fisher.binom.two.sided(rep0, rep1, eps_sub=.01, alpha=0.05)
rep_equiv<-c(1, rep_equiv$odds.ratio, rep_equiv$equiv.CI[1], rep_equiv$equiv.CI[2], "Republican")

third0<-reg_duringAVR$third_party[reg_duringAVR$BSEP25==0]
third1<-reg_duringAVR$third_party[reg_duringAVR$BSEP25==1]
third_equiv<-fisher.binom.two.sided(third0, third1, eps_sub=.01, alpha=0.05)
third_equiv<-c(1, third_equiv$odds.ratio, third_equiv$equiv.CI[1], third_equiv$equiv.CI[2], "Third Party")

white0<-reg_duringAVR$white[reg_duringAVR$BSEP25==0]
white1<-reg_duringAVR$white[reg_duringAVR$BSEP25==1]
white_equiv<-fisher.binom.two.sided(white0, white1, eps_sub=.01, alpha=0.05)
white_equiv<-c(1, white_equiv$odds.ratio, white_equiv$equiv.CI[1], white_equiv$equiv.CI[2], "White")

black0<-reg_duringAVR$black[reg_duringAVR$BSEP25==0]
black1<-reg_duringAVR$black[reg_duringAVR$BSEP25==1]
black_equiv<-fisher.binom.two.sided(black0, black1, eps_sub=.01, alpha=0.05)
black_equiv<-c(1, black_equiv$odds.ratio, black_equiv$equiv.CI[1], black_equiv$equiv.CI[2], "Black")

hispanic0<-reg_duringAVR$hispanic[reg_duringAVR$BSEP25==0]
hispanic1<-reg_duringAVR$hispanic[reg_duringAVR$BSEP25==1]
hispanic_equiv<-fisher.binom.two.sided(hispanic0, hispanic1, eps_sub=.01, alpha=0.05)
hispanic_equiv<-c(1, hispanic_equiv$odds.ratio, hispanic_equiv$equiv.CI[1], hispanic_equiv$equiv.CI[2], "Hispanic")

asian0<-reg_duringAVR$asian[reg_duringAVR$BSEP25==0]
asian1<-reg_duringAVR$asian[reg_duringAVR$BSEP25==1]
asian_equiv<-fisher.binom.two.sided(asian0, asian1, eps_sub=.01, alpha=0.05)
asian_equiv<-c(1, asian_equiv$odds.ratio, asian_equiv$equiv.CI[1], asian_equiv$equiv.CI[2], "Asian")

age18230<-reg_duringAVR$age1823[reg_duringAVR$BSEP25==0]
age18231<-reg_duringAVR$age1823[reg_duringAVR$BSEP25==1]
age1823_equiv<-fisher.binom.two.sided(age18230, age18231, eps_sub=.01, alpha=0.05)
age1823_equiv<-c("OR 2018", age1823_equiv$odds.ratio, age1823_equiv$equiv.CI[1], age1823_equiv$equiv.CI[2], "Age 18-23")


age24300<-reg_duringAVR$age2430[reg_duringAVR$BSEP25==0]
age24301<-reg_duringAVR$age2430[reg_duringAVR$BSEP25==1]
age2430_equiv<-fisher.binom.two.sided(age24300, age24301, eps_sub=.01, alpha=0.05)
age2430_equiv<-c("OR 2018", age2430_equiv$odds.ratio, age2430_equiv$equiv.CI[1], age2430_equiv$equiv.CI[2], "Age 24-30")

age31400<-reg_duringAVR$age3140[reg_duringAVR$BSEP25==0]
age31401<-reg_duringAVR$age3140[reg_duringAVR$BSEP25==1]
age3140_equiv<-fisher.binom.two.sided(age31400, age31401, eps_sub=.01, alpha=0.05)
age3140_equiv<-c("OR 2018", age3140_equiv$odds.ratio, age3140_equiv$equiv.CI[1], age3140_equiv$equiv.CI[2], "Age 31-40")

age41600<-reg_duringAVR$age4160[reg_duringAVR$BSEP25==0]
age41601<-reg_duringAVR$age4160[reg_duringAVR$BSEP25==1]
age4160_equiv<-fisher.binom.two.sided(age41600, age41601, eps_sub=.01, alpha=0.05)
age4160_equiv<-c("OR 2018", age4160_equiv$odds.ratio, age4160_equiv$equiv.CI[1], age4160_equiv$equiv.CI[2], "Age 41-60")

age61plus0<-reg_duringAVR$age61plus[reg_duringAVR$BSEP25==0]
age61plus1<-reg_duringAVR$age61plus[reg_duringAVR$BSEP25==1]
age61plus_equiv<-fisher.binom.two.sided(age61plus0, age61plus1, eps_sub=.01, alpha=0.05)
age61plus_equiv<-c("OR 2018", age61plus_equiv$odds.ratio, age61plus_equiv$equiv.CI[1], age61plus_equiv$equiv.CI[2], "Age 61+")

data_equiv<-rbind(dem_equiv, rep_equiv, third_equiv, white_equiv, black_equiv, hispanic_equiv, asian_equiv, age61plus_equiv, age4160_equiv,  age3140_equiv, age2430_equiv, age1823_equiv)
data_equiv<-data.frame(data_equiv)
data_equiv<-data_equiv %>%
  rename(full=X1, estimate=X2, lower=X3, upper=X4, name=X5)

data_equiv$estimate<-as.numeric(data_equiv$estimate) 
data_equiv$lower<-as.numeric(data_equiv$lower) 
data_equiv$upper<-as.numeric(data_equiv$upper) 
data_equiv$year<-"OR 2018"
data_equiv$full<-"Full Implementation Period"

equivalence_tests<-data_equiv


###Sample tests 
dem0<-sample_OR2018_reg2018$democrat[sample_OR2018_reg2018$BSEP25==0]
dem1<-sample_OR2018_reg2018$democrat[sample_OR2018_reg2018$BSEP25==1]
dem_equiv<-fisher.binom.two.sided(dem0, dem1, eps_sub=.01, alpha=0.05)
dem_equiv<-c(1, dem_equiv$odds.ratio, dem_equiv$equiv.CI[1], dem_equiv$equiv.CI[2], "Democrat")

rep0<-sample_OR2018_reg2018$republican[sample_OR2018_reg2018$BSEP25==0]
rep1<-sample_OR2018_reg2018$republican[sample_OR2018_reg2018$BSEP25==1]
rep_equiv<-fisher.binom.two.sided(rep0, rep1, eps_sub=.01, alpha=0.05)
rep_equiv<-c(1, rep_equiv$odds.ratio, rep_equiv$equiv.CI[1], rep_equiv$equiv.CI[2], "Republican")

third0<-sample_OR2018_reg2018$third_party[sample_OR2018_reg2018$BSEP25==0]
third1<-sample_OR2018_reg2018$third_party[sample_OR2018_reg2018$BSEP25==1]
third_equiv<-fisher.binom.two.sided(third0, third1, eps_sub=.01, alpha=0.05)
third_equiv<-c(1, third_equiv$odds.ratio, third_equiv$equiv.CI[1], third_equiv$equiv.CI[2], "Third Party")

white0<-sample_OR2018_reg2018$white[sample_OR2018_reg2018$BSEP25==0]
white1<-sample_OR2018_reg2018$white[sample_OR2018_reg2018$BSEP25==1]
white_equiv<-fisher.binom.two.sided(white0, white1, eps_sub=.01, alpha=0.05)
white_equiv<-c(1, white_equiv$odds.ratio, white_equiv$equiv.CI[1], white_equiv$equiv.CI[2], "White")

black0<-sample_OR2018_reg2018$black[sample_OR2018_reg2018$BSEP25==0]
black1<-sample_OR2018_reg2018$black[sample_OR2018_reg2018$BSEP25==1]
black_equiv<-fisher.binom.two.sided(black0, black1, eps_sub=.01, alpha=0.05)
black_equiv<-c(1, black_equiv$odds.ratio, black_equiv$equiv.CI[1], black_equiv$equiv.CI[2], "Black")

hispanic0<-sample_OR2018_reg2018$hispanic[sample_OR2018_reg2018$BSEP25==0]
hispanic1<-sample_OR2018_reg2018$hispanic[sample_OR2018_reg2018$BSEP25==1]
hispanic_equiv<-fisher.binom.two.sided(hispanic0, hispanic1, eps_sub=.01, alpha=0.05)
hispanic_equiv<-c(1, hispanic_equiv$odds.ratio, hispanic_equiv$equiv.CI[1], hispanic_equiv$equiv.CI[2], "Hispanic")

asian0<-sample_OR2018_reg2018$asian[sample_OR2018_reg2018$BSEP25==0]
asian1<-sample_OR2018_reg2018$asian[sample_OR2018_reg2018$BSEP25==1]
asian_equiv<-fisher.binom.two.sided(asian0, asian1, eps_sub=.01, alpha=0.05)
asian_equiv<-c(1, asian_equiv$odds.ratio, asian_equiv$equiv.CI[1], asian_equiv$equiv.CI[2], "Asian")

age18230<-sample_OR2018_reg2018$age1823[sample_OR2018_reg2018$BSEP25==0]
age18231<-sample_OR2018_reg2018$age1823[sample_OR2018_reg2018$BSEP25==1]
age1823_equiv<-fisher.binom.two.sided(age18230, age18231, eps_sub=.01, alpha=0.05)
age1823_equiv<-c("OR 2018", age1823_equiv$odds.ratio, age1823_equiv$equiv.CI[1], age1823_equiv$equiv.CI[2], "Age 18-23")


age24300<-sample_OR2018_reg2018$age2430[sample_OR2018_reg2018$BSEP25==0]
age24301<-sample_OR2018_reg2018$age2430[sample_OR2018_reg2018$BSEP25==1]
age2430_equiv<-fisher.binom.two.sided(age24300, age24301, eps_sub=.01, alpha=0.05)
age2430_equiv<-c("OR 2018", age2430_equiv$odds.ratio, age2430_equiv$equiv.CI[1], age2430_equiv$equiv.CI[2], "Age 24-30")

age31400<-sample_OR2018_reg2018$age3140[sample_OR2018_reg2018$BSEP25==0]
age31401<-sample_OR2018_reg2018$age3140[sample_OR2018_reg2018$BSEP25==1]
age3140_equiv<-fisher.binom.two.sided(age31400, age31401, eps_sub=.01, alpha=0.05)
age3140_equiv<-c("OR 2018", age3140_equiv$odds.ratio, age3140_equiv$equiv.CI[1], age3140_equiv$equiv.CI[2], "Age 31-40")

age41600<-sample_OR2018_reg2018$age4160[sample_OR2018_reg2018$BSEP25==0]
age41601<-sample_OR2018_reg2018$age4160[sample_OR2018_reg2018$BSEP25==1]
age4160_equiv<-fisher.binom.two.sided(age41600, age41601, eps_sub=.01, alpha=0.05)
age4160_equiv<-c("OR 2018", age4160_equiv$odds.ratio, age4160_equiv$equiv.CI[1], age4160_equiv$equiv.CI[2], "Age 41-60")

age61plus0<-sample_OR2018_reg2018$age61plus[sample_OR2018_reg2018$BSEP25==0]
age61plus1<-sample_OR2018_reg2018$age61plus[sample_OR2018_reg2018$BSEP25==1]
age61plus_equiv<-fisher.binom.two.sided(age61plus0, age61plus1, eps_sub=.01, alpha=0.05)
age61plus_equiv<-c("OR 2018", age61plus_equiv$odds.ratio, age61plus_equiv$equiv.CI[1], age61plus_equiv$equiv.CI[2], "Age 61+")

data_equiv<-rbind(third_equiv,  rep_equiv, dem_equiv, asian_equiv, black_equiv, hispanic_equiv, white_equiv,  age61plus_equiv, age4160_equiv,  age3140_equiv, age2430_equiv, age1823_equiv)
data_equiv<-data.frame(data_equiv)
data_equiv<-data_equiv %>%
  rename(full=X1, estimate=X2, lower=X3, upper=X4, name=X5)

data_equiv$estimate<-as.numeric(data_equiv$estimate) 
data_equiv$lower<-as.numeric(data_equiv$lower) 
data_equiv$upper<-as.numeric(data_equiv$upper) 
data_equiv$year<-"OR 2018"
data_equiv$full<-"Birthday Window"


equivalence_tests<-rbind(equivalence_tests, data_equiv)
equivalence_tests$name <- factor(equivalence_tests$name, levels=unique(equivalence_tests$name))

saveRDS(equivalence_tests, file="Oregon/Object Files 2018/OR18_equiv.rds")


#ggplot(equivalence_tests, aes(x = estimate, y = name)) +
#  labs(x = "Equivalence range", y = "") +
#  geom_segment(aes(x = lower, xend = upper, yend = name),
#               size = 4, color = "grey90") +
#  geom_point() +
#  xlim(0.9, 1.1) +
#  theme_light()+
#  facet_grid(full ~ year)

sink()


